\documentclass[Orbiter Technical Reference.tex]{subfiles}
\begin{document}

\section{Distributed Vessel Mass}
\label{sec:dist_vessel_mass}
\textbf{Martin Schweiger}\\
\textbf{September 27, 2005}


\subsection{Introduction}
A point mass $m$ placed at position $\vec{r}_0$ in a gravitational field $\vec{g}(\vec{r})$ experiences a force $\vec{F}_G = m\vec{g}(\vec{r}_0)$.
For an extended object with a density distribution $\rho(\vec{r})$ the resulting force can be obtained by integrating over its volume $V \subset \mathbb{R}^3$:
\begin{equation*}
\vec{F}_G = \int_V \vec{g}(\vec{r})\rho(\vec{r}) d\vec{r}
\end{equation*}
For numerical calculations it is sometimes useful to discretise the object into a rigid system of point masses $m_i$ whose relative positions are defined by their barycentric coordinates $\vec{s}_i$. Then,
\begin{equation*}
\vec{F}_G = \sum_i m_i \vec{g}(\vec{s}_i+\vec{r}_\text{CG})
\end{equation*}
where $\vec{r}_\text{CG}$ is the position of the barycentre. For the calculation of the linear force $\vec{F}_G$ Orbiter makes the assumption $\vec{g}(\vec{s}_i + \vec{r}_\text{CG}) = \vec{g}(\vec{r}_\text{CG})$, i.e. the gravitational field is homogeneous over the volume of the object. This approximation is justified when calculating the gravitational force on a spacecraft which is small compared to its orbital radius vector, $|\vec{s}_i| \ll |\vec{r}_\text{CG}|$. With this assumption, we arrive back at the expression for a point mass:
\begin{equation*}
\vec{F}_G = \vec{g}(\vec{r}_\text{CG}) \sum_i m_i = m \vec{g}(\vec{r}_\text{CG})
\end{equation*}
However, an inhomogeneous potential will also induce an angular moment $\vec{M}_G$ in an extended object, and this can generally not be neglected.
In the continuous case, $\vec{M}_G$ is given by
\begin{equation}\label{eq:cont_torque}
\vec{M}_G = \int_V \vec{g}(\vec{r}) \rho(\vec{r}) \times (\vec{r}-\vec{r}_\text{CG}) d\vec{r},
\end{equation}
and after discretisation this becomes
\begin{equation}\label{eq:torque}
\vec{M}_G = \sum_i m_i \vec{g}(\vec{s}_i + \vec{r}_\text{CG}) \times \vec{s}_i
\end{equation}

\subsection{Symmetric gravitational potential}\label{sec:sympot}
If we assume that
\begin{equation}\label{eq:gravpot}
\vec{g}(\vec{r}) = -GMr^{-3}\vec{r},
\end{equation}
i.e. the central body is a sphere of mass $M$ with homogeneous density distribution, and further that $|\vec{r}| = |\vec{r}_\text{CG} + \vec{s}| \gg |\vec{s}|$ over the volume $V$ of the spacecraft, then we can approximate~\cite{wertz1978}:
\begin{equation}\label{eq:radapprox}
r^{-3} = (\vec{r}\cdot\vec{r})^{-3/2} =
\left\{ r_\text{CG}^2 \left[ 1 + \frac{2\vec{r}_\text{CG}\cdot\vec{s}}{r_\text{CG}^2} + \frac{\vec{s}^2}{r_\text{CG}^2} \right] \right\}^{-3/2} \approx
r_\text{CG}^{-3} \left[ 1 - \frac{3 \vec{r}_\text{CG}\cdot\vec{s}}{r_\text{CG}^2} \right]
\end{equation}
Substituting Eqns.~\ref{eq:gravpot} and \ref{eq:radapprox} into Eq.~\ref{eq:cont_torque} leads to
\begin{equation}
\vec{M}_G = \frac{3GM}{r_\text{CG}^3} \int_V (\hat{r}_\text{CG} \times \vec{s}) (\hat{r}_\text{CG} \cdot \vec{s}) \rho(\vec{r}) d\vec{r}
\end{equation}
If the vectors $\vec{s}$ and $\hat{r}_\text{CG}$ are expressed in the vessel reference system, then $\vec{M}_G$ can be written as
\begin{equation}
\vec{M}_G = \frac{3GM}{r_\text{CG}^3} [(\mat{L}\hat{r}_\text{CG}) \times \hat{r}_\text{CG}]
\end{equation}
where $\mat{L}$ is the vessel's inertia tensor expressed in the same frame. (Note that Orbiter currently assumes that the vessel frame of reference is orientated so that $\mat{L}$ is diagonal.)


\subsection{Discrete point mass systems}
Orbiter implements gravity gradient torque as discussed in Section~\ref{sec:sympot}, assuming that the vessel's inertia tensor is known. In this section we give an alternative method that describes the vessel as a rigid system of point masses. This method is not currently implemented in Orbiter.

\subsubsection{2-point systems}
In the simplest case, a vessel can be expressed as a rigid system of two point masses. This allows to simulate an angular moment as a result of a gradient in the gravitational potential $\vec{g}(\vec{r})$.

\setlength{\unitlength}{1mm}
\begin{picture}(100,50)(0,35)
\put(40,70){\circle*{3}}
\put(40,70){\vector(0,-1){25}}
\put(45,58){\makebox(0,0){$\vec{g}_\text{CG}$}}
\put(42,74){\makebox(0,0){CG}}
\put(40,70){\vector(-3,1){20}}
\put(30,77){\makebox(0,0){$\vec{s}_1$}}
\put(40,70){\vector(3,-1){41}}
\put(62,66){\makebox(0,0){$\vec{s}_2$}}
\put(19,77){\vector(0,-1){30}}
\put(15,62){\makebox(0,0){$\vec{g}_1$}}
\put(82,56){\vector(0,-1){15}}
\put(86,48){\makebox(0,0){$\vec{g}_2$}}
\put(19,77){\circle*{2}}
\put(14,77){\makebox(0,0){$m_1$}}
\put(82,56){\circle*{2}}
\put(87,56){\makebox(0,0){$m_2$}}
\end{picture}

\noindent We want to calculate the angular moment induced by the difference of the gravitational fields at $\vec{s}_1$ and $\vec{s}_2$. Let the ratio of masses be denoted by $m_1/m_2 = \gamma$. Then $\vec{s}_2 = -\gamma \vec{s}_1$, and the forces acting on the two mass points are
\begin{equation*}
\vec{F}_1 = m_1 \vec{g}_1,
\qquad \vec{F}_2 = m_2 \vec{g}_2 = \frac{m_1}{\gamma} \vec{g}_2.
\end{equation*}
The angular moment is obtained by adding both components,
\begin{equation}\label{eq:moment}
\begin{split}
\vec{M}_G &= \vec{F}_1 \times \vec{s}_1 + \vec{F}_2 \times \vec{s}_2 \\
&= m_1 \vec{g}_1 \times \vec{s}_1 - \frac{m_1}{\gamma} \vec{g}_2 \times \gamma \vec{s}_1 \\
&= m_1 \left(\vec{g}_1 - \vec{g}_2\right) \times \vec{s}_1
\end{split}
\end{equation}
which is now expressed as a function of the local field gradient $\vec{g}_1 - \vec{g}_2$.

\subsubsection{Numerical implementation}
The field difference is usually small compared to the magnitude of the field acting on the vessel:
\begin{equation*}
|\vec{g}_1 - \vec{g}_2| \ll |\vec{g}_{1,2}|
\end{equation*}
resulting in a significant loss of precision when the field difference in Eq.~\ref{eq:moment} is calculated directly.

To avoid this problem, we can simplify Eq.~\ref{eq:moment} if the field $\vec{g}(\vec{r})$ can be condsidered to be generated by a single point mass $M$ at position $\vec{R}_0$ relative to CG.
Let $\vR{1} = \vR{0} - \vec{s}_1$ and $\vR{2} = \vR{0} - \vec{s}_2$ be the position of $M$ relative to $m_1$ and $m_2$, and let $\vec{s} = \vec{s}_2-\vec{s}_1$.

\setlength{\unitlength}{0.5mm}
\begin{picture}(100,120)(-50,-35)
\put(34,62){\makebox(0,0){CG}}
\put(40,66.5){\circle*{3}}
\put(19,77){\circle*{2}}
\put(81,46){\circle*{2}}
\put(19,77){\vector(2,-1){61}}
\put(19,77){\vector(0,-1){100}}
\put(40,67){\vector(0,-1){89.5}}
\put(81,46){\vector(0,-1){69}}
\put(51,66){\makebox(0,0){$\vec{s}$}}
\put(20,-28){\makebox(0,0){$\vR{1}$}}
\put(41,-28){\makebox(0,0){$\vR{0}$}}
\put(81,-28){\makebox(0,0){$\vR{2}$}}
\put(13,77){\makebox(0,0){$m_1$}}
\put(88,46){\makebox(0,0){$m_2$}}
\put(19,46){\line(1,0){61}}
\put(15,60){\makebox(0,0){$h$}}
\put(19,50){\line(1,0){4}}
\put(23,46){\line(0,1){4}}
\put(19,46){\makebox(4,4){$\cdot$}}
\end{picture}

\noindent Since $\nR{0} \gg |\vec{s}|$, we can find the difference $h = \nR{1} - \nR{2}$ of the distances between the point masses $m_1$ and $m_2$ from $M$ by
\begin{equation*}
h = |\vec{s}| \cos \sphericalangle (\vec{s},\vR{0})
= \frac{\vec{s} \vR{0}}{\nR{0}}
\end{equation*}
We can now write the field difference $\vec{g}_1 - \vec{g}_2$ as
\begin{equation*}
\vec{g}_1 - \vec{g}_2 =
GM \left( \frac{\vR{1}}{\nR{1}^3} - \frac{\vR{2}}{\nR{2}^3} \right) 
\end{equation*}
and by substituting $\vR{2} = \vR{1} - \vec{s}$ and $\nR{2} = \nR{1} - h$, and omitting higher-order terms,
\begin{equation*}
\begin{split}
\frac{\vec{g}_1 - \vec{g}_2}{GM} &= \frac{\vR{1}}{\nR{1}^3} - \frac{\vR{1} - \vec{s}}{(\nR{1}-h)^3} \\
&= \frac{h (- 3\nR{1}^2 + 3\nR{1} h - h^2) \vR{1} + \nR{1}^3 \vec{s}}
{\nR{1}^3 (\nR{1}-h)^3} \\
&= \frac{3h(h - \nR{1})\vR{1} + \nR{1}^2\vec{s}}
{\nR{1}^4 (\nR{1} - 3h)} + O(2)
\end{split}
\end{equation*}
Substituting into Eq.~\ref{eq:moment} and utilising $\vec{s} \parallel \vec{s}_1$ leads to
\begin{equation*}
\vec{M}_G \approx 3 G M m_1 \frac{h(h-\nR{1})}{\nR{1}^4(\nR{1}-3h)}
\vR{1} \times \vec{s}_1
\end{equation*}

\subsubsection{Multi-point systems}
For objects composed of more than two mass points, the formulation must be somewhat extended. We split the gravitational potential acting on mass point $m_i$ into a barycentric and a perturbation component: $\vec{g}(\vec{r}_\text{CG} + \vec{s}_i) = \vec{g}_\text{CG} + \vec{\gamma}_i$.
Then Eq.~\ref{eq:torque} becomes
\begin{equation}\label{eq:multitorque}
\begin{split}
\vec{M}_G &= \sum_i (\vec{g}_\text{CG} + \vec{\gamma}_i) \times m_i \vec{s}_i\\
&= \vec{g}_\text{CG} \times \sum_i m_i \vec{s}_i + \sum_i \vec{\gamma}_i \times m_i \vec{s}_i\\
&= \sum_i \vec{\gamma}_i \times m_i \vec{s}_i
\end{split}
\end{equation}
since the definition of the barycentre demands $\sum_i m_i \vec{s}_i = 0$.

The perturbation components are calculated in the same way as in the 2-point problem, assuming that the gravitational potential is generated by a single point mass $M$ at position $\vec{R}_0$ in barycentric coordinates of the spacecraft, with $|\vec{R}_0| \gg |\vec{s}_i|\;\forall i$.
Given $\vec{R}_i = \vec{R}_0 - \vec{s}_i$, the difference $h_i$ between the distances of $m_i$ and the CG from $M$ is given by
\begin{equation*}
h_i = |\vec{s}_i| \cos\sphericalangle (-\vec{s}_i, \vec{R}_0)
= -\frac{\vec{s}_i \vec{R}_0}{|\vec{R}_0|}
\end{equation*}
Then as before,
\begin{equation*}
\begin{split}
\frac{\vec{\gamma}_i}{GM} &= \frac{\vec{R}_0 - \vec{s}_i}{(|\vec{R}_0|+h_i)^3} - \frac{\vec{R}_0}{|\vec{R}_0|^3} \\
&= -\frac{h_i(3|\vec{R}_0|^2 + 3|\vec{R}_0| h_i + h_i^2) \vec{R}_0 + |\vec{R}_0^3| \vec{s}_i}{(|\vec{R}_0|+h_i)^3 |\vec{R}_0|^3} \\
&= -\frac{3h_i(|\vec{R}_0|+h_i) \vec{R}_0 + |\vec{R}_0|^2 \vec{s}_i}{|\vec{R}_0|^4 (|\vec{R}_0| + 3h_i)} + O(2)
\end{split}
\end{equation*}
and inserting into Eq.~\ref{eq:multitorque} leads to
\begin{equation*}
\vec{M}_G = -\frac{3GM}{|\vec{R}_0|^4} \vec{R}_0 \times \sum_i
\frac{m_i h_i(|\vec{R}_0|+h_i)}{|\vec{R}_0| + 3h_i} \vec{s}_i
\end{equation*}

\subsection{Damping}
The model defined above has equilibrium states ($\vec{M}_G = 0$) for the attitudes $\vec{s} \parallel \vR{0}$ (vessel axis aligned with radius vector) and $\vec{s} \perp \vR{0}$ (vessel axis perpendicular to radius vector). Only the first of these is stable because an attitude perturbation will generate a torque in the opposite direction, leading to an undamped oscillation around the equilibrium attitude.

Orbiter allows to introduce a damping term
\begin{equation*}
\vec{M}_D = -\alpha\vec{\omega}_G
\end{equation*}
where $\vec{\omega}_G$ is the angular velocity induced by torque $\vec{M}_G$, and $\alpha$ is a user-defined damping coefficient. The physical source for the damping term may be the deformation of the vessel by tidal forces, or redistribution of liquid propellants.

\end{document}
